import os
from datetime import datetime, timedelta
from skimage import measure
# Evalscript para calcular NDWI
evalscript_ndwi = """
//VERSION=3
function setup() {
return {
input: [{
bands: ["B03", "B08"] // Green & NIR
}],
output: {
bands: 2
}
};
}
function evaluatePixel(sample) {
return [sample.B03, sample.B08];
}
"""
year = 2023
# Lista para almacenar las áreas de cada mes
areas = []
# Crear una figura y ejes para una cuadrícula de 3x4
fig, axs = plt.subplots(3, 4, figsize=(20, 15))
for month in range(1, 13):
start_date = f"{year}-{month:02d}-01"
if month == 12:
end_date = f"{year + 1}-01-01"
else:
end_date = f"{year}-{month + 1:02d}-01"
ndwi_request = SentinelHubRequest(
evalscript=evalscript_ndwi,
input_data=[
SentinelHubRequest.input_data(
data_collection=DataCollection.SENTINEL2_L1C.define_from(
"s2l1c", service_url=config.sh_base_url),
time_interval=(start_date, end_date),
mosaicking_order=MosaickingOrder.LEAST_CC,
)
],
responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
bbox=pardo_bbox,
size=pardo_size,
config=config,
)
ndwi_array = ndwi_request.get_data()
image = ndwi_array[0]
# Extraer las bandas
Green = image[:, :, 0]
NIR = image[:, :, 1]
# Calcular el NDWI
ndwi = (Green - NIR) / (Green + NIR)
# Crear una máscara para delimitar las zonas de agua
mask = ndwi < 1
ndwi_masked = np.where(mask, ndwi, np.nan)
# Etiquetar los objetos en la imagen binarizada
labeled_ndwi, num_features = measure.label(mask, background=0, return_num=True)
# Calcular las propiedades de los objetos
properties = measure.regionprops(labeled_ndwi)
# Definir un tamaño mínimo (esto depende del tamaño de tu embalse y de la resolución de la imagen)
min_size = 1000 # Ajusta este valor según sea necesario
# Crear una máscara para filtrar los objetos por tamaño
large_objects_mask = np.zeros_like(mask, dtype=bool)
for prop in properties:
if prop.area >= min_size:
large_objects_mask[labeled_ndwi == prop.label] = True
# Aplicar la máscara a la imagen NDWI original
filtered_ndwi = np.where(large_objects_mask, ndwi, np.nan)
# Calcular el área total en metros cuadrados
pixel_area = resolution ** 2 # Área de un píxel en metros cuadrados
water_area = (np.sum(large_objects_mask) * pixel_area) / 10000
areas.append(water_area)
# Obtener los ejes correspondientes a la cuadrícula 3x4
ax = axs[(month-1) // 4, (month-1) % 4]
# Mostrar la imagen NDWI filtrada con el colormap personalizado
cax = ax.imshow(filtered_ndwi, cmap="Blues")
# Título
ax.set_title(f'{start_date[:7]}')
# Ajustar el espacio entre los gráficos
plt.tight_layout()
# Mostrar la cuadrícula de imágenes
plt.show()
# Crear la serie temporal del área cubierta por agua
months = [datetime(year, month, 1).strftime("%b") for month in range(1, 13)]
plt.figure(figsize=(10, 6))
plt.plot(months, areas, marker='o', linestyle='-', color='b')
plt.xlabel('Meses')
plt.ylabel('Área cubierta por agua (ha)')
plt.title('Variación del área cubierta por agua en 2023')
plt.grid(True)
# Guardar la gráfica de la serie temporal
plt.savefig("monthly_images/NDWI_area_variation_2023.png")
# Mostrar la gráfica de la serie temporal
plt.show()
